Slow rotation of a superfiuid trapped Fermi gas 
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The moment of inertia, 0, is one of the possible observables for the experimental determination 
whether a trapped Fermi system has reached the BCS transition or not. In this article we investigate 
in detail the temperature dependence of below the critical temperature T c . Special care is taken 
to account for the small size of the system, i.e., for the fact that the trapping frequency hco is of the 
same order of magnitude as the gap A. It is shown that the usual transport approach, corresponding 
to the leading order of an expansion in powers of fi, is not accurate in this case. It turns out that O 
does not change rapidly if T becomes smaller than T c , but it rather decreases slowly. Qualitatively 
this behavior can be explained within the two-fluid model, which again corresponds to the leading 
order in h. Quantitatively we find deviations from the two-fluid model due to the small system size. 

PACS numbers: 67.,03.65.Sq,05.30.Fk 



I. INTRODUCTION 

Since the first observation of Bose-Einstcin condensa- 
tion of magnetically trapped bosonic atoms |l], ||, |J| it has 
become clear that ultra-cold trapped atomic gases pro- 
vide an excellent tool to study quantum effects in systems 
which are almost visible to the naked eye. For example, 
quantized vortices in the Bose condensate were created 
by stirring the Bose condensate with the help of a laser 
beam |Q . Also the quantum pressure related to the Pauli 
principle could be observed in gases of trapped fermionic 
atoms [| , which proves that temperatures well below 
the degeneracy temperature can be reached. 

If it was possible to trap two spin states of a fermionic 
isotope with attractive interaction, and to cool the sys- 
tem below the critical temperature T c , one could study 
the BCS transition to the superfiuid phase. Unlike the 
transition of a Bose gas to the Bose-Einstein condensate, 
the BCS transition of a Fermi gas almost does not change 
the density profile of the atomic cloud 0. However, there 
are other observables which may allow to distinguish be- 
tween the normal-fluid and the superfiuid phase. In a 
preceding paper || the moment of inertia was proposed, 
since it is much smaller in the superfiuid phase than in 
the normal- fluid phase (see also Ref. ||). Another ob- 
servable changing from one phase to the other are the 
frequencies of collective modes ||, [To[ For example, 
the frequency of the so-called "scissors mode" , an oscilla- 
tion of the symmetry axis of the cloud with respect to the 
symmetry axis of the trap, is closely related to the mo- 
ment of inertia || . Recently one more observable for the 
detection of the BCS transition was proposed, namely 
the change of the deformation of the cloud during the 
expansion of the system when the trapping potential is 
switched off @. 

The moment of inertia of a superfiuid gas of trapped 
fermionic atoms at zero temperature was evaluated for 
the first time in Ref. J8| in close analogy to the calculation 
of the moment of inertia of superfiuid nuclei |]l3| . This 
derivation was very similar to the one given by Migdal 
more than 40 years ago Ej], except that everything was 



reformulated in phase space in terms of Wigner trans- 
forms. In the present article we will generalize the calcu- 
lation of Ref. |8| to the case of non-zero temperature. In 
addition, we will give a derivation which further clarifies 
certain points which in Ref. |^| may have been passed 
over rather quickly. 



In addition to the temperature dependence of the mo- 
ment of inertia, we will address an interesting question 
which is relevant already at zero temperature. In nu- 
clear physics it is well known that the moment of inertia 
of superfiuid nuclei is much smaller than the rigid-body 
value, but still higher than the value corresponding to a 
purely irrotational motion, and that the currents in ro- 
tating nuclei have both rotational and irrotational com- 
ponents Jl3| . The same behavior is found in trapped 
Fermi gases at zero temperature ||. In contrast to this, 
the ordinary hydrodynamical or transport equations for 
superfluids at zero temperature, which can be derived 
from the H — > limit of the time-dependent Hartree- 
Fock-Bogoliubov (TDHFB) equation g| [H], 0, || |l| 
allow only for a purely irrotational motion. We will work 
out this difference and discuss the limits of validity of the 
hydrodynamical description. 



The article is organized as follows: In Sect. [n] we give 
a brief review of the formalism, mainly in order to re- 
call some definitions and to clarify our notation. In Sect. 



Ill we derive the expression for the density matrix of 
the slowly rotating system within linear-response theory. 
This is the generalization of the calculation of Ref. ||] to 
non-zero temperatures. In Sect. IV we again derive the 
linear response of the density matrix, but now using the 
leading order of the h expansion of the TDHFB equation. 
In Sect. [v| we show numerical results for the moment of 
inertia obtained within both formalisms as a function of 
temperature and interprete the results and their differ- 
ences. Finally, in Sect. VI, we summarize and draw our 
conclusions. 
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II. BRIEF REVIEW OF THE FORMALISM 

Before considering the rotating superfluid trapped 
Fermi gas, we will briefly review the equilibrium case. 
Our intention is to explain our notation and conventions. 
Detailed discussions of the subject can be found in many 
articles @, |, || and textbooks [fnj. 

In this article we assume for simplicity that equal num- 
bers of atoms with two spin projections a =T,i arc 
trapped in a spin-independent harmonic potential 



Vb(r) = 



(1) 



=xyz 



If the density of the trapped system is very low, the atom- 
atom interaction can be approximated by a zero-range 
interaction with a coupling constant g proportional to 
the s-wave scattering length. Due to the Pauli principle 
only atoms with opposite spin projections can interact in 
this way. Under these assumptions the hamiltonian takes 
the form 



H = Jd 3 r[ £ ^(r)( 



h 2 v 2 

2m 



Vo(r)) Mr) 



^!(r)v4(r)V T (r)^(r)l . (2) 



The mean-field potential corresponding to this inter- 
action reads 

V(r) = V (r) - gp{v, r) = Voir) - gp(r) , (3) 

where we have used the following notation for the non- 
local density matrix: 



p(r,r') = (V4(r')^(r)> = WWiM) 



(4) 



(Note that with this definition the local part of the den- 
sity matrix, p(r) = p(r, r) corresponds to the density per 
spin state.) In the presence of pairing correlations, the 
pairing gap is given by the gap equation 



A(r) = gn(r, r) 
where the pairing tensor has been defined as 

K (r,r')HV>i(r>T(r)>- 



(5) 



(6) 



It will turn out that the self-consistent solution of Eq. 
(||) is divergent as a consequence of the zero-range inter- 
action. In the literature several ways how to regularize 
this divergence can be found |, f| §(J §§, but in fact 
the technical details of the solution of Eq. (||) are not 
important for our purpose. 

In order to write down the Hartree-Fock-Bogoliubov 
(HFB) equations, which relate the density matrix p and 
the pairing tensor k to the potential V and the gap A, it is 
useful to expand all quantities in a basis of single-particle 



wave functions ip n (r), where n represents all quantum 
numbers except spin, i.e., for an arbitrary operator A: 



A nn , = /dVdV<^(rV„,(r')A(r,r'). 



(7) 



Expressing the field operators ipa( r ) and ^J-(r) m terms 



of annihilation and creation operators a„ 
recover the usual definitions 

Pun' (^n't^ 71 ?) ' 



and al„, we 

(8) 
(9) 

The index n' in Eq. (||) denotes the time-reversed state 
characterized by ipn> (r) = ip* n , (r) . We need also the ma- 
trix elements h nn i of the grand-canonical (mean-field) 
single-particle hamiltonian (i.e., of the single-particle 
hamiltonian minus the chemical potential p) 



P 

2m 



+ Vir)-p, 



(10) 



and the matrix elements A nn > of the gap A. For the more 
general case that the hamiltonian is not time-reversal in- 
variant, we introduce the notation A nn i — An'n- If the 
matrices mentioned above are combined as follows: 



K = 
H = 



p —K 
— K* 1 — p 

h A_ 
At -h 



(11) 
(12) 



the HFB equations | 
a 2 x 2 matrix equation, 



can be written in the form of 



[H,Tl} =0. 



(13) 



What is relevant for our purpose is the spectrum of 
the lowest lying quasiparticles, which for a sufficiently 
small gap can be obtained within the BCS approxima- 
tion, which is much simpler than the solution of the full 
HFB equation (|l3|). We choose a basis in which h is 
diagonal, i.e., h nn > — h n S nn '. Then, within the BCS ap- 
proximation, p and k are diagonal, too, and given by 



2E n 



[1 - 2f(E n )] , 



n n = ^[l-2f{E n )] 



(14) 
(15) 
and 



The quasiparticle energies E n — y/h„ + A^ 
the quasiparticle occupation numbers f(E n ) = 
l/[exp(£ , „/r) + 1] are determined by the diagonal 
matrix elements A„ = A n „ alone. If we neglect the 
non-diagonal matrix elements of A, which are irrelevant 
for the excitation spectrum and, apart from that, much 
smaller than the diagonal ones, we can rewrite Eqs. (|l^ ) 
and dl5|) in the compact form 



n = \-^\i-2 m] 



(16) 
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It is evident that the generalized density matrix TZ given 
by Eq. ([l6]) solves the HFB equation (|lj) if h and A are 
assumed to be diagonal. 

For the spherical case (luq x = luq v = ujq z ) and moder- 
ate numbers of particles (N < 10 4 ), the self-consistent 
HFB equation can be solved numerically [^(|. However, 
for the deformed case and large numbers of particles (ex- 
perimentally numbers of the order TV w 10 5 . . . 10 6 have 
been reached), even within the BCS approximation, the 
self-consistent solution becomes numerically intractable. 
Therefore it may be indicated to apply semiclassical ap- 
proximations. Semiclassical methods can become very 
accurate for large numbers of particles, and in addition 
they often allow for a very clear interpretation of the re- 
sults. To that end we will use the Wigner transforms 
of the density matrix p, the pairing tensor k, the single- 
particle hamiltonian h, etc. The Wigner transform of a 
single-particle operator A is defined as 



A(r,p) = J d 3 se' lp ^ h A^ 



s s 
r H — , r 

2' 2 



(17) 



The Wigner transform h(r, p) of the single-particle 
hamiltonian h is particularly simple: It is just the clas- 
sical hamiltonian. We also recall the useful relations 
[At](r,p) = A*(r,p) and [A](r,p) = A(r,-p). One 
advantage of the Wigner transforms in semiclassical cal- 
culations is the product rule for the Wigner transform 
of the product of two operators A and B pi]], directly 
leading to an h expansion: 



[AB] (r, p) = A(t, p) cxp(^) fl(r, p) , (18) 



where the symbol A stands for the Poisson bracket 



A 



i—xyz 



d_d_ 

On Opt 



_d__d_ 

dpi On 



(19) 



From the definition (|17j) it is clear that the local density 
can be written as 



p(r) = p(r,r) 



d 3 p 
(27rft) ; 



P("",p) 



(20) 



As a very simple case we consider the Thomas-Fermi 
(h — ► 0) limit for the density matrix without pairing 
correlations (i.e., A = n — 0) at zero temperature. 
Quantum-mechanically the density matrix is in this case 
just given by the Fermi sea filled up to the Fermi energy 
p, i.e., p = 6(—h). To leading order in h the Wigner 
transform of this expression gives p(r,p) = 6[—h(r,p)]. 
The corresponding (local) density reads 



P(r) 



67T 2 fi 3 

with the local Fermi momentum 



PF {v) = y/2m\n-V{r)]6\n - V(r)} . 



(21) 



(22) 



Eq. ( |2l| ) together with Eq. ([|) can easily be solved self- 
consistently [0, |8|. Since the pairing gaps and tempera- 
tures considered in this article are very small compared 
with the Fermi energy, we will use Eq. (|2l]) also in the 
presence of pairing correlations and at non-zero tempera- 
tures. (The effect of pairing correlations and temperature 
on the density profile p(r) was investigated in Ref. 0). 
Furthermore, in this article we are not interested in the 
details of the density profile p(r). As shown in Ref. [||, 
the self-consistent solution for p(r) can be described to 
good accuracy by approximating the self-consistent po- 
tential V(r) again by a harmonic potential 



V(r) 



V ^ r 2 

i—xyz 



(23) 



with "effective" frequencies u>i > ujQi (we consider an 
attractive interaction, i.e., g > 0) and an appropriately 
readjusted chemical potential p. In the remaining part 
of this article we will use the approximate potential (23]). 

In order to include the pairing correlations, one can 
also use the HFB equation (|l^) in the limit H — * 0: 

[W(r,p),ft(r,p)] = 0. (24) 

This implies that, to leading order in h, at each point r 
the solution TZ(r, p) as a function of p is given by the 
solution for a homogeneous system with the density cor- 
responding to the local density at this point r [Local- 
Density Approximation (LDA)]: 

1 H(r,p) 

2 2£(r,p) 



ft(r,p) 



[l-2f[E(v,p)}), (25) 



with the definition E(r, p) = y/h 2 (r, p) + A 2 (r). In 
terms of the Wigner transform k(t, p), the gap equation 
[Eq. (||)] can be written as 

A W- 9 /(2W K(r ' P) ' (26) 
Inserting the expression for re(r, p) corresponding to Eq. 
( ^5| ) into Eq. (|26|), we obtain the following non-linear 
equation for the gap: 

d 3 p A(r) 



A(r) = g 



-(l-2/[£(r,p)]). (27) 



(27rft) 3 2E(r,p) 

As mentioned before, this equation is divergent and needs 
some regularization (see Refs. ^ ||, ^2| for details). 

Contrary to the Thomas-Fermi approximation for the 
unpaired density matrix, Eq. (pl|), which is valid if the 
potential can be regarded as constant on a length scale 
of the inverse Fermi momentum, the local-density ap- 
proximation in the paired case is valid only if the poten- 
tial is also constant on a length scale of the coherence 
length of the Cooper pairs. This latter condition is often 
not fulfilled. Therefore, in Refs. ||, |25|] an alternative 
semiclassical method for the calculation of the gap has 
been proposed, which, however, results in an average gap 
(more precisely: gap averaged over the Fermi surface) of 
almost the same magnitude as the average gap obtained 
within the local density approximation. 
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III. LINEAR RESPONSE TO A SLOW 
ROTATION 

In this section we will describe the formalism used for 
the calculation of the moment of inertia of a superfluid 
gas of trapped fermionic atoms. Looking at a system 
rotating with angular velocity fl around the z axis, we 
can calculate the moment of inertia from 



e = 



n 



d 3 r d 3 p 
{2irh) 3 



{r x p v - r yPx ) p(r, p) , (28) 



where p(r, p) is the density matrix of the rotating sys- 
tem. Hence the main problem in calculating the moment 
of inertia is to calculate p(r,p), from which also other 
interesting quantities like the current density (per spin 
state) 



jm - / 



d 3 p 
(2tt/i) : 



— P(r,p) 

m 



(29) 



and the velocity field v(r) = j(r)/p(r) can be derived. 

The system is put into rotation by rotating the exter- 
nal trapping potential around the z axis (of course, for 
this purpose the trapping potential must not be axially 
symmetric). In the rotating frame, however, the system 
is still in a stationary state. In this frame, the Hamilto- 
nian receives the additional term 



hi = —hi 



-Q.L? 



(30) 



which for sufficiently small Q can be treated as a pertur- 
bation. This perturbation induces a change of the density 
matrix, p±, and of the pairing tensor, k\. The mean field 
potential is not changed to linear order in il, since L z is 
a time-odd operator. Linearizing Eq. (113), we obtain 



[H,Tli] = - [Hi, 11} 



(31) 



where H and 1Z denote the unperturbed quantities, while 
Hi and IZi refer to the deviations. Assuming that the 
unperturbed quantities p, k, h, and A are diagonal (BCS 
approximation), we can solve Eq. ( |3l"| ) for pi and Ki. 
(This is equivalent to solving the linearized Gorkov equa- 
tions for the normal and anomalous Green's functions at 
equal times; see, e.g., Ref. E3|.) The solution reads: 



Pi nn' 
Kl nn' 



F nn' h li 



pP A A 

r nn' 5 



nn' 



(32) 
(33) 



where (with the short-hand notation p = p„, p' = p n >, 



h = h n , h' = h ri 



etc.) 



pph 

nn' 


_ (p- 


p')(h + 


h') - (k-k')(A + 


A') 


(34) 






E 2 - E' 2 




J 


F pA _ 

nn' 


(p + p 


-1)(A 


+ A') + (k 


+ K')(h 


+ h') 


(35) 






E 2 - E' 2 






nn' 


(p- 


P')(A- 


A') + (/s- 


K')(h- 


h') 


(36) 






E 2 - E' 2 






pnA _ 
nn 


(i-P 


- P ')(h 




-k')(A 


-A') 


(37) 






E 2 - E 12 







In practice, Eq. (|33j) is an integral equation, since the 
change of the gap, Ai, on the r.h.s. is related to the 
change of the pairing tensor, m, by the gap equation. In 
analogy to Eq. (^6|) the gap equation for the perturbed 
quantities reads 



Ai(r) 



d 3 p 
(2nhy 



«i(r,p) . 



(38) 



The solution of this integral equation contains some 
subtleties. For example, the divergence appearing in Eq. 
( |38| ) as a consequence of the zero-range interaction has 
to be regularized in the same way as the correponding 
divergence of the unperturbed gap equation (j|^) (see ap- 
pendix), and the derivations of Eq. (4.34) in Ref. ||, or 
the second equation after Eq. (16) in Ref. H] are not very 
explicit about this point. However, these problems can 
be circumvented in the following way |l5|, 18[ jlSfl : Sup- 
pose all single-particle wave functions are multiplied by 
the same local phase exp[ic/)(r)]. Then the HFB equation 
( |l3| ) can be rewritten in terms of the gauge-transformed 
matrices 



H = e^He 



-i'I> 



where 



$ = 



4> o 

-<j) 



(39) 
(40) 



(41) 



We will consider (j> as small, i.e., of the order of the per- 
turbation. Then, to linear order in the perturbation, the 
gauge transformed HFB equation reads 



where 



[H,TZi] = -[Hi,n] , 

TZi = Ki +i[$,K] , 
Hi = Hi+i[$,H}. 



(42) 

(43) 
(44) 



In the latter expression one has to take into account that 
h does not commute with <j). Explicitly, for a hamiltonian 
h of the form ( |Io| ) and a local gap A(r) one obtains 



Ira 



hi = -hi = -QL Z - — (p • [V0(r)] + [V0(r)] • p) , 

(45) 

Ai(r) = Ai(r) + 2*0(r)A(r). (46) 
Together with the gauge-transformed gap equation 

d 3 p 



Ai(r) 



{2nhy 



Ki(r,p). 



(47) 



Eq. (Q) is again a system of integral equations which 
for an arbitrary function <f>{v) is completely equivalent 
to the original one, Eqs. ([5l]) and (|3^). However, since 
the perturbation hi is time-odd, the change of the gap, 
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Ai, is purely imaginary and therefore can be eliminated 
by an appropriately chosen function (f>(r). Physically, 
this choice of 0(r) corresponds to a transformation into 
the local rest frame of the Cooper pairs [Q. In this 
particular gauge the linearized HFB equation reduces to 



Pi nn' 
^1 nn' 

and instead of Eq. ( ff7| ) 
termines the phase 0(r) 



= 9 



r nn' "l i 



(48) 
(49) 



we have an equation which de- 



d 3 p 
(2irhy 



Ki(r,p) . 



(50) 



The 



We now proceed to the evaluation of Eq. ( ]f 
unperturbed quantities p and k entering in F£ n , [Eq. (34)] 



can be rewritten in terms of h and k according to the BCS 
relations Jl4| ) and (|l5j). In addition, as in Ref. ||, we 
replace A„ by its average value at the Fermi surface, A, 
because F ph and all other relevant quantities are strongly 
peaked at ef- This allows us to write F£ n , as a function 
of two energies £ = h n and £' = h n i : 



[l-f(E)-f(E')](A 2 +ti>-EE') 



lf(E) 



2EE'(E + E') 
f(E')}(A 2 +^' 



EE') 



2EE'(E - E') 



(51) 



where 



we have introduced 
A 2 and E' 



the 



abbreviations E — 
In contrast to Ref. 



J? 2 + A 2 and E' = y/C 2 + A 2 . 

B, we will not drop the thermal quasiparticle occupa- 
tion numbers f(E) and f(E'). As described in detail 
in Ref. |8|, the Wigner transform of an expression like 



Eq. ( 48 ) can be evaluated semiclassically in the following 
way. First we rewrite Eq. (Eg) as an operator equation: 



Pi 



d£ d? Fr h (Z, ?)6(h - 0~hiS(h - C) • (52) 



Then we use the Fourier representation for the 5 func- 
tions, i.e., S(h — £) = fdt/(2Trh) exp[(h — £)t/h], and 
obtain 

x e ihT ^ 2h h 1 (t)e ihT / 2H , (53) 
where we have introduced the notation 



ft x (t) = e iht/%i xe -m/n _ 



(54) 



To leading order in h the Wigner transform of the product 
of the three operators in the second line of Eq. ( |53| ) can 
be expressed as the product of their Wigner transforms 
[see Eq. (|l|)]. Then the integral over T gives a S function 
of the form S[h(r, p) — £] and the integral over £ becomes 
trivial. 



However, for the operator product in hi(t) [Eq. J54|)] 
we will not use the product rule. In this sense we re- 
sum certain h corrections to all orders. One can also 
say that, since the Wigner transform of Eq. (HJ) involves 
the classical trajectories (see below), the long-time in- 
formation is preserved. On the other hand, developing 
the Wigner transform of Eq. (Q) with the product rule 
(|l9| ) into powers of h would lead to the Wigner-Kirkwood 
H expansion, which is only valid in the short-time limit 
(see Ref. The different treatment of the operator 

products in Eqs. d53] ) and (|54|) is necessary for the fol- 
lowing reason: The operator hi connects states with an 
energy difference of the order huj. This is small compared 
with the Fermi energy, which is the relevant scale for the 
variable £ [since the result pi(r, p) will be used in inte- 
grals over p], but not necessarily small compared with 
the gap A, which is the relevant scale for the variable e 
[this point will become clearer when we investigate the 
function F^(£ + e/2, £ - e/2) explicitly]. 

In the case of the effective harmonic oscillator potential 
( p3| ) the Wigner transform of Eq. (54) can be calculated 
exactly. The result reads 



[h 1 (t)}(r,p) = h 1 [v cl (r,p;t),p cl (r,p;t)} 



(55) 



where r ci (r, p; t) and p ci (r, p; t) are the classical orbits in 
the potential (E3^ corresponding to the initial conditions 
r ci (r, p; 0) = r and p ci (r, p; 0) = p, which are given by 

rf (r, p; t) = n cos^i) + —5- sin(wii) , (56) 



mw. 



pf (r, p; t) = Pi cos(cJii) - mujin sin(wji) 
Putting everything together, we obtain 



(57) 



p 1 (r,p) = JdeFP h (h(v,p)+ £ -,h{v,p)- 

x j^e- ist / n hi[v cl {v,p-,t),p c \v,p-t)). (58) 

Now we proceed to the calculation of the response of 
Pi to the external perturbation hi, neglecting for the mo- 
ment the reaction of the pairing field to the rotation, i.e., 
the p-'Vcf) terms in Eq. J45|). In Ref. || this contribution 
was called the "Inglis-Belyaev term" p{ B . In this case 
the Fourier transform in the second line of Eq. ( |58| ) [with 
hi replaced by hi = —Cl(r x py — r y p x )] can easily be eval- 
uated with the aid of Eqs. ( p6| ) and J57]). Inserting the 
result into Eq. (E3) and observing that F(£,£') is sym- 
metric under the exchange of its arguments we obtain 
[the arguments of /i(r, p) will be suppressed for brevity] 



Pi (r,p) 



nPy 



r yP* 



hujA 



hiijjL 



+ fr x p y r y p x 



huj- 



,h- 



(59) 
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with the definition 



My i UJ X 



(60) 



To simplify the expression (j59|) further we note that the 
distribution function p(r, p) is changed only in the vicin- 
ity of the Fermi surface, provided the Fermi energy is 
large compared with huj±, A, and T. Formally this can 
be inferred from the fact that F ph ([ + e/2, £ - e/2) as a 
function of £ is strongly peaked at £ = 0, which leads us 
to the approximation 



F' h (£+ e /2,£-e/2) 



m , (6i) 



with 



G(x) = 1 



[d£F ph (t + xA,Z-xA). 



(62) 



At zero temperature the integral in Eq. ( |6^ ) can be evalu- 
ated analytically whereas the terms containing the quasi- 
particle occupation numbers f(E) and /(£") have to be 
integrated numerically. After some manipulations the 
function G(x) can be written as 



G(x) = 



arsinh(ir) 
x\/l + x 2 



with E± = ± xA) 2 + A 2 . Within the approxima- 
tion ( |6"T| ) the change of the density matrix corresponding 
to the Inglis-Belyaev term finally takes the compact form 



p{ B (r,p)=m{h(r,p)} [r xPy (l 
- r v p x (l- 

with 



LU+G- + LU-G+ 

cj + G_ — cj_G 

U) + — LU- 



(64) 



(65) 



Now we will consider also the change of the pairing 
field A, i.e., the phase </>(r). As mentioned before, this 
phase will be determined by Eq. (|50|), where Ki(r, p) is 
obtained from the Wigner transform of Eq. (^9|) . Again 
we replace A„ and A„/ entering in F%£, by the average 
value A, which allows us to express F££, as a function of 
two energies: 



[l-f(E)-f(E')]A(£-?) 
2EE'(E + E') 
lf(E)-f(E>)]A(Z-e) 
2EE'(E-E') 



(66) 



Then the Wigner transform of Eq. (^) can be calculated 
semiclassically as given by Eq. (|5^) with pi and F ph re- 
placed by Ki and F Kh , respectively. As it was the case for 



FP h (€ + e/2,Z-e/2), the function F Kh {£ + e/2, £ - e/2) 
is strongly peaked at £ = 0, and we approximate it by 



F» fc (f+e/2,f-e/2) 



with 



G(x) = — d£F Kn ti + xA,€-xA) 
x ' 



(67) 



(68) 



It turns out that the definitions (|62|) and (|68| ) indeed 
define the same function G(x), which is explicitly given 
by Eq. (|63|). Inserting the Wigner transform of Eq. j49| ) 
into Eq. (|50|), we obtain 



d 3 p 



S[h(r,p)] 



{2irhf 

x J^e- i - t / n h l [v c \v,p ] t),p cl {v,p-t)}. (69) 

To solve this equation for the phase 4>{r) we make the 
ansatz 101 



(r) 



(70) 



Then the second line of Eq. (|69| ) is just the Fourier trans- 
form of h x = -tt(r x p y - r y p x ) - a(r x py + r y p x ), which is 
readily evaluated with the aid of Eqs. ^5^) and (p7|). Due 
to the 5 functions the remaining integrals are trivial, and 
Eq. (KM) finally becomes 







igm 2 p F (r)r x r y 



8ir 2 H 2 A 

x [fiw+u>_(G + + GL) + a(^G+ + w 2 ^)] , ( 7i ) 
which has the solution 

w+w_(G+ + G_) 



= -fi- 



LU 2 G+ +L0 2 _G- 



(72) 



Using this expression we can also calculate the change 
of the original pairing field, Ai: Since the change of 
the gauge-transformed pairing field [Eq. (fl6|)] is zero, the 
original pairing field is modified according to 



Ai(r) = -2iA0(r) = - 



2iA 



amr x r y 



(73) 



Having calculated the phase 4>(r), we can now evaluate 
Eq. ( p8| ) with the full h\, i.e., including in addition to 
the Inglis-Belyaev term [Eq. ((59|)] also the response of 
the density matrix to the p • terms. This second 
contribution to p\ , which we will call p^ 1 , is obtained in 
the same way as discussed above for the first one, and 
the result reads 



pf^T,p) = ad[h(v,p)] 
+ r y p x ( 



r x p y (1 - 
oj+G 



1 



oj+G+ + w_G_ 
o_-w_G_' 



UJ+ 



(74) 
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However, we are not interested in the change of the 
gauge-transformed density matrix, pi , but of the original 
density matrix, p\. According to Eq. ( [43| ) the relation 
between p\ and p\ is given by 



P1 = h-A^p\ = pr + P r i +pr- (75) 

Due to the simple r dependence of </>, the Wigner trans- 
form of the commutator [0, p] is identical to the Poisson 
bracket of the Wigner transforms of 4> and p, i.e. 



pf 2 (r,p) = fi0(r)Ap(r,p) 



d 

d P y 



' dp x 



P(r,p). 



(76) 



As we did before, we will assume that A and T are much 
smaller than the Fermi energy. Therefore we can write 
p(r,p) ~ 0[—h(r,p)] and we obtain 



2 0", P) = -a(r xPy + r y p x )S[h(r, p) 



(77) 



The total effect of the phase </>, i.e., of the reaction of 
the pairing field, on the density matrix, which in Ref. ||] 
was called the "Migdal term" pf 1 , is the sum of the two 



contributions p^ 11 and p\ 



Mi 



Pi (r,P) = -ao[h{r,p)\[r x py ■ 



TyPx ) ■ (78) 



Together with the Inglis-Belyaev term, Eq. (64), and the 
explicit expression for a, Eq. (fT2|), our final result for the 
change of the density matrix reads 



(79) 



pi(r, p) = QS[h(r, p)] r x p y - r y p. 



AG+G-(uj 2 x r x p y - u:lr y p x ) 
ujIG++lg 2 _G- 



Given the change of the density matrix, we can imme- 
diately calculate the current density j(r) [Eq. (p9|)] and 
the velocity field v(r): 



j(r) = (o(r)v(r) = Clp(r)(r x e y - r y e x 

AG + G-{L0 2 x r x e y - uj 2 r y e x 



G. 



(80) 



It is interesting to check explicitly that this current ful- 
fils the continuity equation. In the rotating frame the 
continuity equation reads 



V • j(r) + p(r) - Q(e z x r) • Vp(r) = , 



(81) 



where p(r) = in our case of a stationary rotation. Tak- 
ing the divergence of Eq. ( jSp] ) , we get from the second line 
a contribution proportional to [Vp(r)]'[e z xVV(r)]. This 
is zero, since the gradient of the density in Thomas-Fermi 



approximation [Eq. (|2T|)], Vp(r), is parallel to VV(r). 
Thus, the divergence of the current is equal to the diver- 
gence of the first line of Eq. (|8(]), which exactly fulfils 
Eq. (|si|). Note that the contribution of the Migdal term 
is crucial in order to satisfy the continuity equation. The 
easiest way to see this is to consider the limit A — > oo. 
In this limit we have G± — > 1 and p{ B (r,p) — > 0, which 
implies j (r) 



-> 0. Hence, with the Inglis-Belyaev con- 
tribution alone, Eq. (p 
As observed in Ref. 



cannot be satisfied. 
, the velocity field v(r) describes 
a mixture of rotational motion, corresponding to a veloc- 
ity field proportional to e z x r, and irrotational motion, 
corresponding to a velocity field proportional to V(r x r y ). 
The ordinary rigid rotation is realized if G+ = G_ = 0. 
This is the case if the temperature approaches the criti- 
cal temperature T c , where the gap vanishes [the temper- 
ature dependence of the function G(x) will be discussed 
in Sect. M), but it can also happen at zero temperature 
if A <§; hw±, as discussed in Ref. ||. Purely irrotational 
motion, as it is expected in homogeneous superfluids, is 
reached if G+ = G_ = 1. This is only possible if the 
temperature is very low and if A 3> Hu>±. 

For completeness let us also discuss the change of the 
pairing tensor, «i(r, p), which can be obtained in a way 
completely analogous to the calculation of the change of 
the density matrix, pi(r, p). The result reads 



Ki(r,p) 



2mn 

mA 



UJ+UJ. 



- G . G- 



<G+ 



G- 



-p x p y S[h(r,p) 



- 2ic 



(r)«(r,p). (82) 

Since the last term is of the order h~ x [see Eq. ((to|)] , it has 
been argued that a semiclassical description is possible 
only in the particular gauge where A + Ai is real and 
where this term vanishes 1181, 19(1 . 



IV. SUPERFLUID ROTATION IN TRANSPORT 
THEORY 

The transport or hydrodynamical equations for super- 
fluid systems can be derived by taking the h — > limit of 
the time-dependent Hartree-Fock-Bogoliubov (TDHFB) 
equation @ |lj || 

ihn=[H,n], (83) 

i.e., by replacing the Wigner transforms of the commuta- 
tors in Eq. ( 83 ) by Poisson brackets of the Wigner trans- 
forms [Q |16j, |17|, [l8| , [TfJ . Due to the transformation into 
the rotating frame, we are dealing with a static problem, 
where the TDHFB equation (p3| ) reduces to the HFB 
equation (jl^) . Again we make use of the gauge transfor- 
mation and retain only terms of linear order in the per- 
turbation. Then, if 0(r) is chosen such that Ai vanishes 
[Eq. (^FJ)], the leading order in H of Eq. ( ^2[ ) becomes 



ihhApi + 2Aki = —ihhiAp, 
ihAKpi — 2hk\ — ihhiAn . 



(84) 
(85) 
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In this equation and in the remaining part of this sec- 
tion, h, p, k, etc. denote the Wigner transforms of the 
corresponding operators; the arguments r and p are sup- 
pressed for brevity. 

Let us first study the zero-temperature limit, T = 0. 
In this case the unperturbed quantities are given by p = 
(1- h/E)/2 and k = A/2E [see Eq. @ in the limit 

T — * 0], and it is easy to show that {h\h.p)h = {hKn)A. 
Thus, for A ^ 0, the solution of Eqs. (ph and Q reads 



Pi =0, 



(86) 
(87) 



As we will see, the relation p± = implies that the ve- 
locity field is completely irrotational independent of the 
magnitude of A, which is a well-known property of ho- 
mogeneous superfluid systems at T = 0. 

Now we are going to determine the phase <f>. To that 
end we insert Eq. (^7|) into Eq. (^0|). If we make again 
the ansatz (|70|), we obtain the following equation: 







igh 
2A 



d 3 P (, n . s dp 



(CI - a)r y 



dp \ 
dr x ) 



[Note that in this equation p still refers to the Wigner 
transform of the non-local density matrix, p(r,p).] It 
is clear that in general this equation does not have a 
solution for all r, since the ansatz ( |70[ ) is not general 
enough. But under certain assumptions it turns out that 
this ansatz is sufficient. Firstly, we assume that the gap 
A(r) is either replaced by a constant corresponding to its 
average value at the Fermi surface (as it was done in the 
previous section), or that A(r) is calculated within the 
LDA. In these both cases the function A(r) can formally 
be written as A[y(r)]. Using this, we define the following 
short-hand notation: 

dp _ dp dh dp dA _ A 2 hA dA 

dV ~ dh dV + dA dV ~ ~2E3 + dV ' ^ ' 

which allows us to write Vp = (dp/dV) VV. Secondly, 
as in the previous section, we assume that the potential 
V(r) is a harmonic oscillator. Then Eq. (88) becomes 



ighmr x r y 
2A 



d 3 p dp 



(2nh) 3 dV 
x [SHwl-uD +a(w2+w»)], (90) 



As in the previous section, the phase cj) implies a change 
of the density matrix, pi , due to the inverse gauge trans- 
formation, which to leading order in h reads 



Pi — pi + h(f>Kp . 



(92) 



As we have seen, the first term vanishes. Thus, to linear 
order in the perturbation, Eq. ( |92"|) can be rewritten in 
the following, more suggestive way: 



P(r, P) + Pi (r, p) = p[r, p + fiV^(r)] . 



(93) 



From this equation it follows immediately that the veloc- 
ity field is given by 



v(r) 



-V0(r), 



(94) 



which is completely irrotational. Note that this result 
does not depend on the form of c/)(r) and the approxi- 
mations made to calculate 4>(r). It also does not at all 
depend on the magnitude of A, as long as A ^ 0. It is 
rather a direct consequence of the vanishing of p± , which 
in turn follows immediately from the h — > limit of the 
linearized HFB equations for time-odd perturbations and 
zero temperature. However, as we have seen in the pre- 
vious section, in a small system where hui is of the same 
order of magnitude as A, the velocity field is not irrota- 
tional. Our conclusion is that one should be careful when 
applying transport theory to such systems. 

So far we have considered only the zero-temperature 
limit. In the remaining part of this section we are going 
to consider also the case T > 0. In this case it is difficult 
to solve the coupled Eqs. ( p4| ) and (|85|). However, if we 
in analogy to the previous section assume that the the 
unperturbed gap A is constant, we find the following 
solution for pi and Ki: 



_ (dp Adny df(E) t 



\dh h dh) 



dE 



Kl 



ih~ dn ~ Y 



(95) 
(96) 



If we again make the ansatz ( |70| ) and insert Eq. (96) into 
Eq. ( p0[ ) , we find a — ao as in the zero-temperature case 
[see Eq. (pl|)]. This could have been anticipated from 
the h — > limit of Eq. (|7^), which does not depend on 
the actual value of G(0). Finally we are now going to 
calculate p\. To that end we insert Eqs. ( |95| ) and ( |70| ) 
with a = ao hito Eq. (|92|), and we obtain 



with the solution 



a 



(91) 



Not surprisingly, this result is identical to the h — > 
limit of Eq. ([72]), since for T = we have G(0) = 1 and 
consequently lim^o G± = 1. 



Pi 



-fl ^_ (r x p y 



dE 



a °(^-S (r ^ +r ^ ) - (97) 



Since df(E)/dE and dp/dh are both strongly peaked at 
the Fermi surface, we can make the same approximation 
as in the previous section, i.e., we replace the strongly 
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peaked functions by 6 functions with the appropriate 
strength. Noting that 



lim Gix) = 1 + / 

x^O J 



we can write the result as 



df(E) 
dE ' 
(98) 



Pi = n[l - G(0)]8(h)(r x p y - r y p x ) 

- a Q G(0)S(h)(r xPy + r yPx ) , (99) 

which is in perfect agreement with the h —> limit of 
Eqs. (fj) and ©. 



V. RESULTS AND DISCUSSION 

Using the results for change of the non-local density 
matrix pi(r, p) given in the previous sections, we can now 
calculate the moment of inertia. It should be remembered 
that an ideal Fermi gas at zero temperature behaves like a 
rigid body, i.e., the velocity field is given by v(r) = fle z x 
r. Since the critical temperature for the BCS transition 
is very low, 9 will approach the rigid-body value Qrigid 
for T — > T c . Using the Thomas- Fermi density profile ( |2l| ) 
with the effective harmonic oscillator potential (^), we 
can immediately calculate Q r i g id- The result reads 



4r 1 



rigid 



I2h 3 



UJz-UJy-UJz 



(100) 



In terms of Q r igid the moment of inertia of the superfluid 
system as obtained from pi(r, p) can be written as 



e = e, 



rigid 



8u%J*G+G. 



{ul + ul){ulG++^_G-) 



(101) 



I n th e h — > (transport) 
(|l0l|) reduces to 



limit, where G± -> G(0), Eq. 



9 = 9 



rigid 



l - 



G(0) + G(0)(| 



(102) 



In fact, this formula can be understood very easily. 
The moment of inertia corresponding to the purely irro- 
tational velocity field as it is expected for a large super- 
fluid system at zero temperature, v(r) = — oio^(r x r v ), is 
given by 



9, 



= 9 



( L 

riqid \ 



(103) 



Within the two-fluid model a homogeneous system of 
density p is described as a mixture of a superfluid com- 
ponent of density p s and a normal-fluid component of 
density p n , with p s + p n = p. At T = one has p s = p 
and p n = 0, whereas at T > T c one has p s = and 



p n = p. If this model was correct also for finite systems, 
one would expect that the moment of inertia is given by 



v ' vJ rigid 

P 



yJtrrot • 

p 



(104) 



This would be exactly Eq. (102), if we could identify G(0) 
with /O s /p. In fact, the microscopic calculation of p s for 
a homogeneous system gives [p7| 



1 



Ps= P 



6n 2 mh 3 



df(E) \ 
dE J 



(105) 



with E = v /(p 2 /2m - p) 2 + A 2 . Noting that the in- 
tegrand is peaked at p — pf and remembering p — 
p F /6TT 2 h 3 , we rewrite this as 



df(E) 
dE 



(106) 



with ^ = p /2m — p. As noted in Sect. IV, the r.h.s. of 



this equation is identical to lim^^o G(x), so that we are 
left with 



Ps_ 

p 



G(0). 



(107) 



The previous paragraph can be summarized in the 
statement that the transport approach, corresponding 
to the leading order of the H expansion, reproduces the 
two-fluid model for homogeneous systems. It does not 
give any finite-size corrections, as can be seen from the 
fact that the result does not depend on the trapping fre- 
quencies, except for the purely geometrical dependence 
contained in O r igid and Qirrot- In contrast to this, the 
method described in Sect. |l| is capable to describe the 
different behavior of the system depending on whether 
the trapping frequencies (multiplied by H) are small or 
large compared with the gap A. This dependence is gov- 
erned by the G± factors appearing in Eq. (101), result- 
ing from the long-time behavior of the operator h\ (t) [see 
discussion after Eq. (^3)]. In order to reproduce this be- 
havior within the h expansion, one would have to resum a 
certain class of corrections proportional to huj±/A to all 
orders, in particular if one wants to cover the whole range 
of possible parameters from ftio± < A to hui± A. 

Let us now proceed to a quantitative analysis. In order 
to calculate the moment of inertia 9 as a function of 
temperature, we need the temperature dependence of the 
gap A. As in Ref. J^], we will assume that it is described 
by the same universal function relating A/Aq to T/T c in 
homogeneous matter, where A denotes the gap at T = 
and T c — 0.567Ao This universal function is given by 
the solution of the non-linear equation u& 



f(E) 
E 



(108) 



For completeness it is displayed in Fig. [I]. 

For the calculation of the moment of inertia we also 
need the function G(e/2A), which depends on T via the 



10 



o 




1 I — . — I — - 

A = 10 nK 
0.8 -A = 20nK 
. h -> limit 
"§> 0.6 - ff> x /co v = 0.5 
© 



0.4 



0.2 



o) x /co y = 0.8 



0.2 0.4 0.6 0.8 
T/T„ 



FIG. 1: Universal function for the temperature dependence 
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FIG. 2: Behavior of the function G(e/2A) for different tem- 
peratures between and T c . 



temperature dependence of A discussed above, and via 
the explicit temperature dependence of the function G(x) 
due to the thermal quasiparticle occupation numbers as 
given by Eq. (|63| ) . If only the temperature dependence of 
A was included, G(e/2A) as a function of e would become 
very strongly peaked at e — for T — > T c . However, due 
to the explicit temperature dependence of the function 
G(x), the peak is suppressed and as a function of e the 
function G(e/2A) even becomes more and more flat with 
increasing temperature, as shown in Fig. ||. The decrease 
of G(0) when T approaches T c reflects the decrease of the 
superfluid fraction in the two-fluid model. 

Next we have to specify the parameters of the sys- 
tem. For our comparison we consider, as in Ref. [g, 
583000 6 Li atoms (i.e., 286500 atoms per spin state) in 
a harmonic oscillator potential with average frequency 
hui = h^xUjyLUz) 1 / 3 — 8.21 nK. The corresponding chem- 
ical potential is /x = 983 nK. In order to simulate the 
effect of the self-consistent mean-field potential, the fre- 
quency has been chosen slightly higher than the fre- 
quency of the external trapping potential (huo — 6.9 nK) 



FIG. 3: Moment of inertia as a function of T/T c for small 
(uo x /uo y — 0.8) and large (uj x /ujy = 0.5) deformations in the 
xy plane and two values of Ao (10 and 20 nK). The short- 
dashed lines correspond to the ft — * limit. 



In the experiments the traps are generally very elon- 
gated, i.e. we have a strong deformation a = lo z /lo±, 
where u± — ^jLo x uj y is the average frequencey in the xy 
plane. In our examples we choose a = 1/8. This results 
in a rather high value for the average transverse frequency 
of Huj^ = faj/a 1 / 3 = 16.42 nK. In order to rotate the 
system around the z axis, at least a small deformation 
in the xy plane is necessary, which we parametrize by 
6 = uj x /u;y. (In practice, the rotating deformation of the 
potential can be generated by a laser beam 0].) 

The main uncertainty comes from the gap at zero tem- 
perature, Aq. Note that the coupling constant g does not 
appear explicitly. The moment of inertia depends on the 
interaction only via A, which can be written as a func- 
tion of T and Ao- The value of the critical temperature 
T c = 0.567Ao is still under investigation. In addition, 
the s-wave scattering length a of the atoms, and conse- 
quently g, Ao, and T c , can be tuned in the experiments 
by a magnetic field due to the presence of Feshbach res- 
onances. Therefore we will treat Ao as a free param- 
eter. As a rough estimate, using the scattering length 
a = — 2160eto, where ao is the Bohr radius, one obtains 
that the gap Ao averaged over the Fermi surface is of the 
order of magnitude of 15 nK ||, i.e., of the same order of 
magnitude as the transverse trapping frequency lu±. 

In Fig. U we display the moment of inertia as a func- 
tion of the temperature for two different deformations 
5. The lower curves correspond to a very small deforma- 
tion, 8 = 0.8. In this case the moment of inertia at T = 
is very small. When T approaches T c , the normal- fluid 
component becomes more and more important and con- 
sequently the moment of inertia increases until it finally 
reaches the rigid-body value at T = T c . Qualitatively the 
behavior is similar in the case of a strong deformation in 
the xy plane (upper curves), except that in this case the 
whole curve is shifted upwards, mainly due to the much 
larger value of Qi rro t- The difference between the three 
curves shown for each deformation will be discussed be- 
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FIG. 4: Current distributions j(r x , r y , 0) in the a;y plane (ar- 
bitrary units) for oj x /lo v — 0.8, Ao = 20 nK, and four different 
temperatures: T/T c = 0, 0.4, 0.6, and 1. 



low. 

In order to illustrate the origin of the temperature de- 
pendence of 0, we show in Fig. the current distributions 
for the case S = 0.8 for four temperatures between T = 
and T = T c . One can clearly see the continuous transi- 
tion from the irrotational motion at T = 0, resulting in a 
small angular momentum and therefore a small moment 
of inertia, to the rigid motion at T = T c . 

Now we are going to discuss the differences between 
the three curves shown in Fig. ^ for each deformation. 
The short-dashed lines correspon d to t he results obtained 
within the h — > approach, Eq. ( |102[ ) . The l ong-dashed 
and solid lines were obtained from Eq. ( |l0l| ), i.e., they 
take into account the difference between G± and G(0), re- 
sulting from the long-time behavior of the operator h\{t) 1 
Eq. (p4|). From the definition (65) it is clear that this 
difference is less important for large values of A, and in- 
deed the long-dashed lines, corresponding to A = 20 nK, 
are closer to the h — > results than the solid lines, cor- 
responding to Ao = 10 nK. More precisely, the crite- 
rion for the validity of the h — > approach seems to be 
fko±_ <C Ao rather than fku± -C A, as one might expect. 
This surprising fact can be understood by looking at Fig. 
|2|: Whatever is the actual value of the temperature T [i.e., 
of A(T)], the value of G(e/2A) can always be replaced 
by G(0) if e/2A < 1. 

To show more clearly the non-trivial dependence of 
on A, we show in Fig. [5] the moment of inertia for zero 
temperature as a function of A for the same deformations 
as in Fig. ||. The irrotational limit, indicated by the 
dashed lines, is reached for A — > 00. If A is much smaller 
than fhu- (3.67 nK in the case S = 0.8 and 11.61 nK in the 
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FIG. 5: Moment of inertia for zero temperature as a func- 
tion of A for small (u] x /lu v — 0.8) and large (uo x /u y — 0.5) 
deformations in the xy plane. The dashed lines indicate the 
corresponding irrotational limits. 



case S — 0.5, respectively), the moment of inertia even 
approaches the rigid body value, and the H expansion 
fails completely. For example, in nuclear physics strong 
deviations from the irrotational value are quite common 

Finally let us briefly discuss the question whether the 
moment of inertia is suitable to detect the superfluidity 
in experiments. In principle the moment of inertia can 
be measured directly by measuring the rotational energy 



E ro t — "j^ 2 



(109) 



Since the rotation does not change the potential energy 
(at least not to linear order in 11), the rotational energy 
is equal to the difference of the release energies E re i of 
the rotating system and of the non-rotating system. (The 
release energy E re i is the total energy of all particles after 
the trapping potential has been switched off, i.e., the 
sum of the kinetic energy Eki n and of the interaction 
energy E int of the trapped system.) A disadvantage of 
the direct measurement of E rot is that it requires two 
identical systems, one in rotation and one at rest. As a 
rough estimate we approximate the release energy E re i by 
the kinetic energy E^n of the particles in the "effective" 
harmonic potential (p3|), 



Ekin 2 



d 3 rd 3 p p 2 . 

T2^F^ (r ' p) 



8fl 3 LU x U)yUJ z 



(110) 



Hence, as a function of the average transverse trapping 
frequency luj_ = ^Juj x lo v and the deformation 5 — lo x /lo v 
we obtain 



Er. 



1 + s 2 / n \ 2 

Ekin 3(5 rigid \LU_l_J 



(111) 



Since we used linear response theory, our results arc valid 
only for slow rotations, <C uj±. In particular the angu- 
lar velocity must be small enough in order to avoid the 
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creation of vortices. For an optimistic estimate we choose 
£1 = QAuj±. Since the difference of G between T = T c 
and T = is most pronounced for small deformation (see 
Fig. ||), we choose 5 — 0.8. Using these numbers we find 
Erot/Ekin ~ 0.1 x O /0 rigid, i.e., the moment of inertia 
might indeed be measurable. 



VI. SUMMARY AND CONCLUSIONS 

In this article we have discussed the temperature de- 
pendence of the moment of inertia of a Fermi gas trapped 
in a slowly rotating trapping potential. The assumption 
of a slow rotation allowed us to use linear response theory 
(RPA), but it is clear that in this way certain interest- 
ing effects like the creation of vortices could not be con- 
sidered, since they depend non-linearly on the angular 
velocity fl of the rotation. 

In Sect. [II we derived the density matrix of the ro- 



tating system using a semiclassical method similar to the 
one described in Ref . || , but now taking into account the 
thermal quasiparticle occupation numbers, which were 
neglected in Ref. || and which give rather important 
contributions. One important point is that the method 
takes into account that the energy difference Huj± of the 
states connected by the perturbation hamiltonian (i.e., 
essentially by L z ) is not necessarily negligible in compar- 
ison with the gap A. This leads to a non-trivial behavior 
of the density matrix on hui±/A. These effects can also 
be regarded as finite-size effects, since hcu± vanishes in 
homogeneous systems. 

In Sect. IV we presented an alternative method for the 



calculation of the density matrix, where only the leading 
order of the h expansion is retained. This is equivalent 
to the transport or hydrodynamical approach which is 
often used in the literature || |l0[ |l2). The qualitative 
difference between the results obtained within the two 
app roaches is that the velocity field obtained in Sect. 
[II has irrotational and rotational contributions at all 
temperatures, whereas the transport approach presented 
in Sect. IV gives a purely irrotational velocity field at zero 
temperature, as it is the case in homogeneous systems. 
The dependence on Hlo±/A mentioned above is missed 
within this approach. 

In Sect. [y| we used the density matrices obtained in 
the preceding sections for the calculation of the moment 
of inertia. The result can qualitatively be understood 
within the two-fluid model, which describes the super- 
fluid system as a mixture of a superfluid and a normal- 
fluid component. The density of the normal-fluid com- 
ponent is zero at T = and approaches the total density 
for T — > T c . We have shown that the transport approach 
exactly reproduces this two-fluid model. Somewhat sur- 
prisingly, the condition for the transport approach to be 
valid turns out to be Hoj -C Ao, where Ao is the value 
of the gap at T = 0. This is less restrictive than the 
condition fku <C A, in particular for temperatures near 



Within the transport approach, the moment of inertia 
increases smoothly from the irrotational limit at T = 
to the rigid-body value at T — T c . This is a consequence 
of the increasing density of the normal-fluid component 
of the two-fluid model. If the condition Hoj <C Ao is 
not fulfilled, the behavior is qualitatively similar, but the 
moment of inertia is always larger than it is within the 
transport approach, because in this case the rotational 
contributions to the velocity field are always non-zero 
due to the finite-size effects mentioned above. In both 
cases, the smoothly increasing moment of inertia as a 
function of temperature can be obtained only if the ther- 
mal quasiparticle occupation numbers are properly in- 
cluded in the calculation. It is not sufficient to perform 
a zero-temperature calculation and then replace the gap 
A by the temperature-dependent gap A(T). 

Looking at the size of the error made by neglecting 
the finite-size effects, we conclude that for the trapped 
fermionic atoms, where hu> < Ao, the hydrodynamical 
approach is just at the limit of its applicability. However, 
we would like to point out that there are other physical 
situations, where fko > Ao, and where finite-size correc- 
tions are crucial. For example, the moments of inertia 
of rotating superfluid nuclei (T = 0) have at least twice 
the irrotational value |l3| . Also for the description of su- 
perconducting metallic grains in a weak magnetic field, 
corresponding to a perturbation hi = (e/mc)p • A(r) = 
(e/mc)B z L z (if B is parallel to the z axis) and there- 
fore being formally equivalent to a slow rotation, these 
corrections might be important. 

The method used in Sect. Ill for the semiclassical so- 



lution of the RPA in superfluid systems can also be ex- 
tended to the dynamical case, i.e., to time-dependent 
perturbations. In this way collective excitations of the 
superfluid system, in particular the change of their fre- 
quencies compared with the normal-fluid phase, can be 
described. So far the collective modes in the superfluid 
phase have been studied either within the hydrodynami- 
cal approach |], [l0| or quantum-mechanically for the case 
of spherical symmetry and moderate numbers of particles 
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APPENDIX: ALTERNATIVE DERIVATION OF 
THE MIGDAL TERM 



In Sect. Ill wc derived the change of the pairing field, 
Ai(r), via a gauge transformation. Here we will present 
an alternative method, which is more direct, but also 
somewhat more difficult. We will solve the original inte- 
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gral equation for Ai which one obtains by inserting Eq. 
© into Eq. f|): 



A 1 (r)=AP(r)+Af (r) 



Mi 



d 3 P i ju , ., m 



with 



.IB 



[«f(r >P ) + «f(r,p)]. (A.l) 



— p k " h 

1 nn' — Jr nn'" / 1 nn' > 
pkA a 

r nn' nn' ■ 



nn 



(A.2) 
(A.3) 



The Wigner transforms of these two contributions to Hi 
can be calculated semiclassically as given by Eq. ( |5q) 
with pi replaced by k[ b and nf 1 , respectively, and F^" 
replaced by F' K ' 1 and F kA , resp ectively. 

The first term in Eq. (A.l), A{ B , has already been 
evaluated in Sect. Ill [term proportional to f2 in Eq. (71)] 
with the result 



Af(r) 



igm 2 p F (r)r x r y 
8ir 2 H 2 A 



ftw+w_(G+ + G_). (A.4) 



Now we turn to the evaluation of the second term, A^. 
The explicit expression for -F kA (£,£') reads 



l-2f(E) 1-2/QEQ 

4E 4E' 
[l-f(E)-f(E'M-C) 2 
4EE'{E + E') 

[f(E)-f(E'M~Q 2 
4EE'(E — E') 



(A.5) 



which in analogy to FP h (£,C) and F K ^ (£,£') can be ap- 
proximated by 



(A.6) 

with E = y^ 2 + A 2 . Using this approximation we get 
[the arguments of the functions h(r, p) and E(r, p) are 
omitted for brevity] 

d 3 P n-2f(E) 



Af(r)= 5 



(2irh) 3 



2E 



Ai(r) 



A e -^ Al [r^(r,p;i)] 



(A.7) 



At this stage the disadvantage of the present method as 
compared with the method used in Sect. [II becomes ob- 
vious, since we encounter a divergent integral over d 3 p, 



whereas in Sect. Ill all expressions were finite. This di- 
vergence is the same one which also appears in the gap 
equation ( p7|) for the unperturbed gap in local-density 
approximation. If we assume that this equation is reg- 
ularized in some way, we can use it to get rid of the 
divergence in Eq. (A.7), and we obtain 



Af(r) = Ax(r)+ 9 



d 3 p 



6(h) 



2A/ V2A 



(2ttK) 3 

f dt e -fet/» Al[r ci ( t)]< (A 
J 2nn 



As we will see, the integral equation (A.l) can be solved 
by the ansat z (|7^ ). With this ansatz the Fourier trans- 
form in Eq. ( |A.S| ) can easily be evaluated and we obtain 



Af(r) = A x (r) + 



igm 2 p F (r)r x r y 

87T 2 ft 2 A 



a{uj 2 + G + +uj 2 _G-). 

The coefficient a can now be determined by inserting 
Eqs. flAl| ) and flA~9| ) into Eq. ( [0] ). The solution, of 
course, coincides with Eq. (|72j). 

However, we have to admit that the above arguments 
concerning the divergence in Eq. (|A.7 | ) are a little bit 
hand-waving. For example, Eq. (|27|) (including an ap- 
propriate regularization) is valid only in the local-density 
approximation, and it does not allow for a constant gap 
A, while we have for simplicit y as sumed that A is a con- 
stant in order to derive Eq. (A.7). Such inconsistencies 
do not appear within the formalism presented in Sect. 

ili. 



It remains to show that the Migdal term, calculated as 
the second term of Eq. (p2[), is consistent with the result 
given in Eq. ([78]) . This can be done with the aid of the 
explicit expression for F pA , which turns out to be 



and the Fourier transform of Ai [r ci (r, p; t)]. 



(A.10) 
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